NASA Contractor Report 198331 
ICASE Report No. 96-30 


//v -6- 


ICASE 


A WAVELET-OPTIMIZED, VERY HIGH ORDER 
ADAPTIVE GRID AND ORDER NUMERICAL METHOD 


Leland Jameson 


NASA Contract No. NAS 1-19480 
May 1996 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, VA 23681-0001 

Operated by Universities Space Research Association 



National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23681-0001 



A WAVELET-OPTIMIZED, VERY HIGH ORDER 
ADAPTIVE GRID and ORDER NUMERICAL METHOD 


Leland Jameson 1 
Mitsubishi Heavy Industries, Ltd. 

Advanced Technology Research Center 
8-1, Sachiura, 1-Chome, Kanazawa-ku, 

Yokohama, 236, Japan 
email: lmj@atrc.mhi.co.jp 

ABSTRACT 

Differencing operators of arbitrarily high order can be constructed by in- 
terpolating a polynomial through a set of data followed by differentiation of 
this polynomial and finally evaluation of the polynomial at the point where a 
derivative approximation is desired. Furthermore, the interpolating polyno- 
mial can be constructed from algebraic, trigonometric, or, perhaps exponen- 
tial polynomials. This paper begins with a comparison of such differencing 
operator construction. Next, the issue of proper grids for high order poly- 
nomials is addressed. Finally, an adaptive numerical method is introduced 
which adapts the numerical grid and the order of the differencing operator 
depending on the data. The numerical grid adaptation is performed on a 
Chebyshev grid. That is, at each level of refinement the grid is a Chebyshev 
grid and this grid is refined locally based on wavelet analysis. 
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1 Introduction 


One can argue that high order numerical methods are appropriate for prob- 
lems in 1) the direct numerical simulation of turbulence, see [17], 2) flows 
with shocks and nonlinear physics, see [8] and 3) flows with smooth propagat- 
ing structures such as those encountered in aeroacoustics. Assertion number 
3) is based on convergence properties of the hp-refinement method in finite 
elements, see [20], [1], [11], in which convergence is very fast for high order 
polynomials as long as the function at hand is smooth. In addition, high 
order methods are more efficient for long time integration of unsteady flow 
problems, see [18]. 

This paper introduces a numerical method which combines very high 
order differencing with a wavelet-based grid and order selection mechanism. 
Here very high order differencing will be schemes of order greater than or 
equal to 8, i.e., perhaps 16, 20, or maybe order 32. Such high orders of 
accuracy can produce solutions which are very close to those produced by 
spectral methods. See [3] for a spectral method on arbitrary grids. 

The numerical method introduced here is named the Wavelet- Optimized 
Finite Difference method 2, or WOFD2. WOFD2 is an extension of WOFD, 
see [15], [9], in which wavelets are used for grid refinement, see [16], for finite 
difference schemes. The new method extends this scheme to very high orders 
and also adapts the order of accuracy depending on the data. 

Let us begin by studying various manners in which high order difference 
operators can be constructed. 


2 Generating Difference Equations 

Given a vector of N numbers / how can we get an approximate value of the 
derivative /' at the i — th point and how good will this approximate_value be. 
Generally speaking, the more elements around the i — th point of / that are 
used to approximate f the better the approximation will be. Common finite 
difference formulas are found by fitting a algebraic polynomial of degree q lo- 
cally around the i — th point of a vector / of evenly-spaced elements to obtain 
difference approximations of accuracy q — 1. This section will generalize this 
concept to find the difference equations of arbitrary accuracy on arbitrary 
grids using algebraic, trigonometric, cosine and exponential polynomials. As 
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special cases, one can obtain all the usual finite difference formulas as well as 
the Fourier collocation and Chebyshev collocation spectral differential ma- 
trices. 

Two methods of generating the differencing coefficients will be introduced. 

The first method explains how to set up a system of equations which will have 
as a solution the differencing coefficients. The second method is the deriva- 
tion of differencing coefficients by interpolation. It is this second method 
which is used throughout the paper for the actual generation of difference 
equations. 

2.1 Setting up a Linear System 

The problem is to find a set of coefficients {r*} which combines the raw data 
in a vector / to provide an approximation to a derivative: 

h=right 

f'( x j) = 2 »■*/(**)• (i) 

k=left 

If we require the above equation to be exact for polynomials, algebraic, trigonometric, 
cosine or exponential, then a linear system of equations can be solved to find 
an appropriate set of differencing coefficients {r*}. Let b(x) denote a funda- 
mental basis element from which a basis can generated by taking powers of 
6(x): b(x) = x, b(x) = e tx , b(x) = cos(x), or b(x) = e x . That is, we require 
that the derivative be exact up to a given order N on the numerical grid. 

The system of equations to be solved for a centered differentiation stencil is 
as follows: 

= J2 r k (b(x j+k )) n , (2) 

k=-L 

^From this equation one can generate a system of equations with N re- 
quirements, that N functions be differentiated exactly, and N degrees of 
freedom, the N differencing coefficients r*. If one is near a boundary, then 
the stenciled is biased. Since this type of system is well-known for algebraic 
polynomials, an example for the less well-known trigonometric polynomials 
will be given. 

Consider a trigonometric polynomial on a 3 point centered stencil. The 
first equation simply requires that the derivative of a constant be zero: 

0 = r_i +r 0 + ri. (3) 
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Note that is the same equation as for algebraic polynomials since (x)° = 
(e ,x )°. The next two equations come from requiring that the n = 1 mode is 
differentiated exactly: 


ie ix ° = + r 0 e' x ° + r ie ixi . (4) 

One now obtains the two equations from equating the real and imaginary 
parts. These three equations can be solved for the three coefficients r_ x , r 0 , rq . 

Similarly, one can find the coefficients for higher order schemes by requir- 
ing that more modes be differentiated exactly. Note that no restrictions were 
placed on the grid. Differencing formulas can be found on arbitrary grids as 
easily as they can be found on uniform grids. Also, note that the Fourier 
spectral differentiation matrix can be found from the above procedure by 
requiring that the grid be uniform and that the differencing formulas have 
maximum accuracy on a given grid. That is, if one is working on a grid 
of size 33 then require that the first 16 modes and the zero-th mode are 
differentiated exactly. 

2.2 Interpolation 

A second approach, and the one used in this paper, is to generate differencing 
coefficients by first interpolating a polynomial through a set of data, followed 
by differentiation of this polynomial and evaluated at a grid point. 

The main reason that differentiation was studied with a variety of types of 
differentiation operators was to find out if there was any advantage to using, 
say, trigonometric polynomials to differentiate as opposed to algebraic poly- 
nomials when the function to be differentiated was for example a Gaussian 
pulse. It seemed like an appropriate study to undertake given the current 
research activity in the area or aeroacoustics where one is often confronted 
with the need to computationally propagate some type of wave motion. The 
thought was that perhaps trigonometric polynomials might have some ad- 
vantage at propagating wave motion over the more common algebraic poly- 
nomials. One of the conclusions of this section is that there is no advantage 
and that one should simply use algebraic polynomials for the generation of 
differencing equations. In fact, the only important issues involved with ob- 
taining approximate derivatives is the order of the finite difference operator 
and the density of the numerical grid. 
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The most important reference for this section is [7]. The following four 
subsections will cite the interpolation formulas for the four types of interpo- 
lation, and hence differentiation, considered in this section. 

2.2.1 Algebraic Polynomials 

Interpolation with algebraic polynomials is probably the most common form 
of interpolation, and it is from this type of interpolation that common uni- 
form grid finite difference methods can be found. Using the following for- 
mula one can find the finite difference coefficients for an arbitrary grid and 
of arbitrary order. One simply fits the polynomial to the data, followed by 
differentiation of the polynomial, and finally one evaluates the polynomial 
at the point of interest. The well-known Lagrange interpolation formula for 
algebraic interpolation is, 


n n 

A A x ) = n ( x - x k)t n ( x ; - **)• ( 5 ) 

kzzO,k^j k=O t k^j 

Aj(x k) = 6jk For given values w 0 , u>i, ..., u>„, the polynomial 

n 

Pn( x ) - WkAk(x). (6) 

Jfc=0 

in P n and takes on these values at the points x,-: 

Pn(xk) = Wk, (7) 


for k = 0, 1, ..., n. 

2.2.2 Trigonometric Polynomials 

As seen from the previous section, one can also generate difference operators 
by using trigonometric functions as the fundamental interpolation elements. 
The following is the appropriate Lagrange- type interpolation formula, see [7]: 
For —7 r < x 0 < xi < ... < x 2n < tt then 

2 n ^ 2n ^ 

r i( a; )= II sin-(x-x*)/ IJ sm-{xj - x k ). (8) 

*=0,*^ Z k=0,k^j Z 
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The function, 


(9) 


2 n 

T(x) = '£w k T k (x) 

k= 0 

is the unique solution of the interpolation problem, 

T(x k ) = w k , (10) 

for k = 0, l,...,2n. Again, one can derive finite difference coefficients by 
interpolating to a function, followed by differentiation of the interpolation 
polynomial and evaluation at the point of interest. The following section will 
prove that such difference equations obey order properties just as the usual 
difference equations derived from algebraic polynomials do. 

2.2.3 Cosine Polynomials 

The comparable Lagrange-style interpolation formula for cosine polynomials 
is the following, see [7]. 

Given n + 1 distinct points 0 < xq < x\ < ... < x n < ft- Set 

n n 

Cj(x) = P (cosx — cosxt)/ (cosx_, — cosxjt). (11) 

k=0,k^j k=O y k^j 

Then Cj is a cosine polynomial of order < n, Cj(x) = a * cos (kx), for 
which Cj(xk) — Sjk. Given n+ 1 distinct values wo , iwi, ..., w n there is a unique 
cosine polynomial of order < n, C(x), for which C(xk) = Wk , k = 0, 1, ..., n. 
It is n 

C(x) = ± Wk Ck(x). (12) 

k = 0 

Note that ^(7(x)|o, w = 0 since ^Cjfe(x)|o,a- = 0 for all k. For this rea- 
son, difference operators based on cosine polynomials will not be considered 
in general, but will be compared in a later section to Chebyshev spectral 
methods. As above, the differencing coefficients are found by first fitting 
the trigonometric polynomial to the data, followed by differentiation of the 
polynomial and finally evaluation at the point of interest. 
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2.2.4 Exponential Polynomials 

The final polynomial to be tested is the exponential polynomial, 


«,•(*)= n («•-«•*)/ n 

k=0,k^j k=0,k^j 

(13) 

where the interpolation polynomial is, 


E(x ) = 53 u>kE k (x), 

k = 0 

(14) 

and 


E(x k ) = w k . 

(15) 


2.3 Truncation Error and Differentiation Accuracy 

The purpose of this section is to illustrate algebraically that one obtains dif- 
ferentiation order-of-accuracy properties for all four types of differentiation 
operators which are similar to the standard order-of-accuracy properties ob- 
tained with the usual algebraic interpolation. In short, if one interpolates 
an N order polynomial then one obtains a reduction in differentiation error 
of (i)^ when the density of the grid is doubled. This order of accuracy is 
obtained regardless of the type of polynomials which are used. 

Recall that the remainder for algebraic polynomial interpolation is, see 

P], 

/w - m*) = ( x - xo)( (;;;jr (l ~ i - > / iw,( o. w 

where £ lies between the smallest and the largest x,. The following sec- 
tion will show that a similar expression can be obtained for any polynomial 
constructed from powers of a given function. 

2.3.1 Truncation Error for Interpolation by Powers 

There are some subtle issues concerning a general proof of truncation error 
and accuracy for interpolation by a polynomial constructed from the powers 
of a general function g(x ), see [4]. The following demonstration will illustrate 
the essentual algebraic steps that one follows to obtain accurcy while avoiding 
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the subtle issues. In short, let the polynomial element g{x) and the function 
to be approximated f(x) be “well-behaved”. 

Let 

p(x) = Y,a k (g(x)) k 


k=0 


be the polynomial which interpolates /(x) at xo,xi, ...,x n , p(x,) = /(x,), 
then, 

/(X) P(X) — A(n4.U/t\ (!') 




where 


4>(x) = (g(x) - g(x 0 ))(ff(x) - <?(xi))-..(ff(x ) - j(i n )), (18) 

and where £ lies between the smallest and the largest x t . 


Demonstration of Truncation Error: Note that much of this demon- 
stration is the same as that which can be found in a standard numerical 
analysis text for the remainder term in algebraic interpolation, see [5]. 
Define H(z ) such that 

H(z) = f(z ) - p(z) - R(x)<f>(z), (19) 


where R(x) is defined such that H(x) = 0. Note that H(xi ) = 0, for i = 
0, ...,n since p(x t ) = /(x,) and <^(x,) = 0. /.From Rolle’s theorem it follows 
that there exists a point £ in the interval defined by the smallest and largest 
x,-’s such that f7^ n+1 ^(£) = 0. This implies, 


R(x) 


< X n+1 >(0 


( 20 ) 


Now put this back into the expression for H(z ) and set z = x to get, 


f(x)-p(x) 


/ (n+1) (0-p (n+1) (0 

<£ (n+1) (0 


^(x). 


This is the desired expression. / /. 


( 21 ) 
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Note that in the above demonstration that if the polynomial is algebraic 
that p( n+1 \z) = 0 and 4>^ n+1 \z) = (n + 1)!, but for a general <?(x) these 
two functions are just a measure of the smoothness of the basic interpolation 
element, g(x). /( n+1 )(£) still remains a measure of the smoothness of the 
function one is interpolating to, and <f>{x ) is a function dependent on the grid 
distribution. 

2.3.2 Differentiation Accuracy 

The primary interest here is to understand the behavior of the derivative 
operators derived from the various types of interpolation outlined above as 
the grid is refined. That is, 

f(x)-p'(x) = Q(0<f>'(x), (22) 

where Q(£) = — - and will dictate the behavior as the 

grid is refined. It will be shown that the behavior of <j>'{x) is essentially 
independent of the basic interpolation element g(x ) and depends only on the 
order of the interpolation. 

Demonstration of Accuracy: 

Let h denote the smallest spacing in the numerical grid, then 

<f>\x 0 ) = Ch n + h.o.t. (23) 

where n is the highest power in the interpolation polynomial, and the point 
xo is an arbitrary grid point inside the interpolation stencil. 

Demonstration: First of all, 

4>'{x o) = g'(x 0 )(g(x 0 ) - ^(xjX^xo) - ^(x 2 ))...( 5 (x 2; ) - g(x n )). (24) 

Expand about zero the function p(x), 

g(x) = g( 0) + g r (0)x + /'(0)x 2 /2 + ... (25) 

and examine the difference <7(x 0 ) — ^(xi) 

g(*o) -g{xi) = g'(0){x o - Xi) +g"(0)/2(xl - x$)+g"'(0)/6(xl - x?) + ... (26) 
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Without loss of generality let one of the points be zero, say x 0 = 0, to get, 
s( 0) - j(x,) = S '(0)(-x,) + S "(0)/2(-x?) + s'"(0)/6(-x?) + .... (27) 


or, 

5 (0) - g(x 0 = -x 1 ( y £ (28) 

m= 1 

and one can see that the first term in the difference is linear. If Xq is not zero 
then one obtains the factorization, 

x>-p» = (x-y)(5>y'- 1 - i ), (29) 

1=0 


and hence, 


00 ( 0 ) m_1 

sM - SM = (*« - x,)( E E (30) 

m=l m - k = 0 

It should be clear that the first term in the difference 5 ( 2 : 0 ) — 5 ( 2 : 1 ) is the 
linear term and that doubling the grid such that between every two points 
another point is placed is, therefore, halves the distance to a first order 
approximation. For each of the differences x t — xj there exists a constant c,j 
such that the difference, 

x^ Xj — Cijh (81) 

can be expressed in terms of the smallest grid spacing h. Let C = YVj=i coj 
then it is, therefore, clear that 

<f>’{x 0 ) = Ch n + h.o.t. (32) 


//• 

Consider the following special cases which includes all the polynomial 
types discussed above: 

• g{x ) — x, 5 ^ m ^( 0 ) = 0, for m / 1 

• g(x) = e x , 5 (m) (0) = 1, Vm 
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g(x) = e**, 5 ( m )(0) = i m 


• g(x ) = cos(x), 0< m )(O) = 0, for m odd and ^”*^(0) = (— l) m / 2 , for m 
even 

For a simple illustration, consider the following two examples of algebraic 
and exponential interpolation. 

Algebraic Interpolation 

Consider the simple case of interpolating an algebraic quadratic polyno- 
mial p 2 (x) to a function f(x ) at the grid points x 0 < X\ < x 2 : p 2 (x;) = /(x t ), 
i = 0,1,2. The remainder term for some £, x 0 < £ < x 2 , is 

/(x) - p 2 (x) = (x - x 0 )(x - Xj)(x - x 2 )^/ (3) (£). (33) 

Now, differentiate and evaluate at x = xi to get, 

f'(x i) - P 2 (x i) = (xj - x 0 )(xi - x 2 )^/ (3) (^) = Ck 2 /&(£), (34) 

where h = x\ — x 0 = x 2 — x\. If the grid is evenly-spaced then the differences 
(x, — Xj) are some integer multiple of the smallest difference which one can 
denote by h. If one doubles the number of grid points then each of the 
distances (x, — Xj) becomes half as large and the accuracy for this quadratic 
example will be 2. 

The General Statement for Algebraic Interpolation 

In general, one can expect that algebraic polynomial interpolation with 
a polynomial of order n, p n (x), will produce a differencing operator of order 
also of order n. This can be seen from the portion of the truncation error 
which depends on the grid distribution: 

n 

— x t ). (35) 

»=o 

This product contains n + 1 terms. After differentiation and evaluation at a 
point Xk the product will contain n terms, 

n— 1 

(36) 

t=0 
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When the grid density is doubled by adding a point between every two points, 
then each distance ( Xk — x,) will decrease by a factor of 2. The product of 
these factors will be a multiple of (|) n and hence the accuracy of n is achieved. 

Exponential Interpolation 

If, on the other hand, p 2 (x) is now an exponential polynomial then after 
differentiation and evaluation one gets, 

/'(* ,) - j 4(*,) = - e-)(e- - e**)i/< 3 >K) (37) 

Without loss of generality, let X\ = 0. Then the product, 

/'( 0) - p' 2 ( 0) = (e*» - e Xo )(e Xl - e x ’) (38) 

becomes 

/'(0) - p^O) = (*o + xl/2 + ...)(x 2 + x\!2 + ...), (39) 

and one can see that if the distance to zero is halved and hence Xo and x 2 are 
divided by 2 that the leading order terms in each of the above parenthesis 
dictates that the error will be reduced by 4 to a first order approximation. 

Therefore, one can expect the accuracy to behave as it does for algebraic 
polynomials. That is, doubling the grid points will decrease the error by 
(1/2) 2 to first order. 

The General Statement for Exponential Interpolation 

In general, one can expect that exponential polynomial interpolation with 
a polynomial of order n, p n (x), will produce a differencing operator of order 
also of order n which is the same result as for algebraic interpolation. As can 
be seen from the above example, differentiation and evaluation at a point x* 
will produce a leading order term which is a product of n terms, 

JI(zfc-a:;), (40) 

t=0 

and accuracy of order n is achieved. 
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Table 1: A Comparison of Errors for Various Types of Differentiation Oper- 
ators 


2.4 A Numerical Check of Accuracy 

This subsection will verify that the differentiation operators which are gener- 
ated from algebraic, trigonometric, and exponential polynomials all exhibit 
the same order property, which depends only on the order of the polynomial 
interpolation. The difference operators will be tested on the function, 

*)« I 

2 + cos(2x) 

defined on [0, tt\. f(x) is chosen because it is periodic but not exactly a 
trigonometric or algebraic polynomial. Table (1) illustrates the order prop- 
erty All of the errors are L 2 - 



A general question arises, is there any advantage to using, say, a trigono- 
metric polynomial for the generation of difference equations over, say, the 
usual algebraic polynomial? From this study the answer appears to be no. 
The essence depends on the ability of the interpolating polynomial to locally 
approximate the function at hand. If the function at hand is not exactly a 
trigonometric or algebraic polynomial, as is likely, then there is no advantage 
for either approach. Such an issue is important when one is considering the 
propagation of, say, a pulse in the application of aeroacoustics. A pulse will 
locally be neither a algebraic or trigonometric polynomial. In short, a wave, 
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or pulse, can not be propagated accurately if it can not first be differentiated 
accurately, and it can not be differentiated accurately if it can not first be 
approximated accurately. 


3 High Order Methods 

Spectral collocation methods are often given the probable misnomer of “in- 
finitely accurate”. In a manner consistent with finite difference methods, 
the accuracy of spectral collocation methods will be assigned the accuracy 
of AT - 1 when applied on a grid of N points. 

This section will begin by connecting spectral collocation methods to fi- 
nite difference methods. That is, spectral methods will be viewed from the 
point of view of the maximum finite difference method on a given grid. Fol- 
lowing the comments on this connection, a case will be made for applying very 
high order algebraically generated finite difference operators on Cheybshev 
grids or, equivalently, applying very high order cosine polynomial generated 
finite difference operators on a uniform grid. This second process of using 
cosine polynomials will require a mapping of the independent variable from, 
say, x to cos(x), but is exactly equal to applying algebraic polynomials on 
Chebyshev grids. 

3.1 Spectral Collocation = Maximum Order Finite 
Difference 

On a numerical grid of N points one can fit a polynomial with N degrees-of- 
freedom through all of the data. If this polynomial is algebraic and if the grid 
distribution is Chebyshev, x, = cos(^), then one can build the Chebyshev 
collocation differentiation matrix. On the other hand, if the polynomial is 
trigonometric and if the grid is uniformly distributed then one can build 
the Fourier spectral collocation differentiation matrix. One can, therefore, 
define in the physical space a spectral method to be a method which uses the 
maximum size polynomial for approximation and differentiation for a given 
grid size. 

Another way to see this is, suppose one has a numerical grid of 16 points 
and a 4-th order difference operator on a 5 point stencil. Now reduce the 
number of grid points to 8. The difference operator is still 4-th order. Now 
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Table 2: Fourier Spectral Collocation Applied to sin’s 


reduce the number of grid points to 5. The difference operator is still 4- 
th order accurate. This is spectral accuracy. That is, spectral accuracy of 
collocation methods on finite grids is N — 1 where N is the number of grid 
points. 

3.1.1 A Numerical Check 

Let us consider the above statements numerically. A Fourier collocation spec- 
tral method on a grid of 9 points is a differencing mechanism with exactly 
9 degrees-of-freedom, and hence, is designed to differentiate exactly 9 func- 
tions exactly: 1, cos(fcx), and sin (kx), for k = 1,2, 3, 4. Table (2) is meant 
to illustrate two points: i) the maximum frequency which is differentiated 
exactly is 4.0, and ii) the poor performance on non-integer frequencies. 

Likewise, a Chebyshev collocation spectral method on a grid of 9 points 
is designed to differentiate 9 functions exactly: x k for k = 0, ...,8. Table (3) 
is designed to illustrate the same two points at Table (2). The table begins 
with the function x 5 . 

Note that if the function being differentiated is not exactly an integer 
frequency, e lkx or x k , then, say for Chebyshev, differentiating x 5 - 5 is no bet- 
ter or worse than the result for differentiating sin(5.5x). The point that is 
trying to be made is that the differentiation accuracy of spectral collocation 
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Grid 

Pts 

Poly 

Order 

l 2 

Error 

9 

x 5 

2.2910 -25 

9 

x 5 - 5 

9.2610“ 4 

9 

X 6 

6.4910 -25 

9 

X 6 - 5 

2.7710 -3 

9 

x 7 

2.3910 -24 

9 

X 7 - 5 

1.7910 -i 

9 

X 8 

6.9610- 24 

9 

X 8 - 5 

4.2910 -1 

9 

X 9 

2.8310° 


Table 3: Chebyshev Spectral Collocation Applied to Polynomials 


methods on a finite grid of size N is accurate with the accuracy of N — 1, 
and one can not expect that a wave-like pulse will be transmitted better 
with a Fourier spectral method than with a Chebyshev spectral. The only 
important issue is the dimension of the space and the boundary conditions: 
use Chebyshev for non-periodic boundary conditions and Fourier for period 
boundary conditions. 

3.2 Very High Order Finite Differencing 

Now suppose that build a series of algebraically generated finite difference 
operators of increasing accuracy and test these difference operators on the 
function sin(2x). The grid size will be fixed at 33 points. The first two lines of 
Table (4) illustrate the effect of the Runge phenomenon, see [5]. The change 
in the error from periodic boundary conditions on a uniform grid to non- 
periodic boundary conditions on a uniform grid is from 10 -27 to 10~ 21 . In 
addition, note that applying an algebraic polynomial with periodic bound- 
ary conditions yields a result comparable to applying a trigonometric, i.e. 
Fourier spectral collocation, polynomial with periodic boundary conditions. 
Furthermore, one can observe the Runge phenomenon with trigonometric 
polynomials just as one observes it with algebraic polynomials when the 
boundary conditions are not periodic. Note that 128 bit arithmetic is being 
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Grid 

Order 

l 2 

Error 

Periodic 

Grid even 

Pts 

of Acc 

Error 

Ratio 

BC’s 

or Cheby 

33 

32 

1.5210“ 27 


yes 

even 

33 

32 

5.0710 -21 


no 

even 

33 

18 

7.4910 -17 


no 

Cheby 

33 

20 

1.1710 -18 

64.0 

no 

Cheby 

33 

22 

1 .7310“ 20 

67.6 

no 

Cheby 

33 

24 

2.4310" 22 

71.2 

no 

Cheby 

33 

26 

3.4910 -24 

69.6 

no 

Cheby 

33 

28 

5.7510 -26 

60.7 

no 

Cheby 

33 

30 

1.3010 -27 

44.2 

no 

Cheby 

33 

32 

00 

OO 

<X> 

i— ‘ 

o 

1 

to 

00 

1.46 

no 

Cheby 


Table 4: Finite Difference Accuracy Approaching Spectral Accuracy 


used. ^From line 3 to the bottom of the table the order of accuracy is in- 
creased from 18 to the maximum, i.e., spectral, accuracy of 32 is obtained. 
When one tests the accuracy of a finite difference operator one doubles the 
grid and sees the error decrease as (^) n where n is the accuracy of the scheme. 
This comes from the truncation error which will produce a factor of the form 
(Ai) n . In Table (4), it is the number n which is being increased while Ax 
remains constant. 

Compare Table (4) to Table (5) in which a Chebyshev collocation method 
on an increasing grid size is tested on sin(2x). Note that in the following table 
that both Ax is decreasing and n is increasing in the expression (Ax) n as 
one proceeds down the table. The final line of Table (5) is the Chebyshev 
method on a grid of 33 points which produces a result comparable the result 
in Table (4) on a grid of 33 points. The numbers are not exactly the same 
because, first of all, all calculations are near machine accuracy, and, second, 
the differencing coefficients are calculated in different ways. 
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Table 5: Chebyshev Spectral Collocation of Increasing Order 
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3.3 Chebyshev Spectral Methods and Cosine Polyno- 
mials 

This subsection will review Chebyshev spectral methods and the equivalency 
with cosine polynomials. Chebyshev approximation can be seen as approxi- 
mation by algebraic polynomials, 

Tq{x) = 1, (42) 

Ti(x) = x, 

Ti{x) = 2x 2 — 1, 

T 3 (x) = 4x 3 — 3x, 
or as approximation by a cosine series, see [7], 

n n 

T n (x ) = cos(n arccos x) = cos(n0) = ^<z,(cos0) ? = ^a ? x 9 , (43) 

9=0 g =0 

for some set {a g } and where x = cos(0). If one now chooses a numerical grid 
defined as Xj = cos(^), j = 0, ...,N then one obtains T n (x 3 ) = cos(^) and 
the pseudospectral Chebyshev method, see [10], [21], and [2]. 

A Chebyshev spectral method involves approximating a f un ction, /(x), 
by interpolating an algebraic polynomial to point values /(x,) where the grid 
points are given by the ‘Chebyshev’ grid points x t - = cos(0,). An equivalent 
process is to interpolate a cosine polynomial to the evenly-spaced point values 
of the angle 9 X and to consider /(x) to be evaluated on the uniform grid of 
the angle variable 0 t , /(x,) becomes /(0,). That is, see [10], if the Chebyshev 
series for f(x) is 

OO 

p f( x ) = 9( x ) = a kTk(x), (44) 

fc = 0 

then the expansion coefficients {a*} can be found in two equivalent ways, 

flfc = / f{ x )Tk(x)(l — x 2 )~ 1/,2 dx = f f(cos6) cos k$d$ (45) 

7TCfc J - 1 7 re* Jo 

By this transformation of the independent variable one can perform Cheby- 
shev spectral methods on a uniform grid or one can build arbitrarily high 
difference operators on a uniform grid which have stability characteristics 
equivalent to the usual Chebyshev spectral method. 
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3.4 High-Order Differencing on Chebyshev Grids 

Chebyshev spectral methods work very well for non-periodic problems pre- 
cisely because the truncation error for a Chebyshev polynomial is equil-ripple. 

As shown above, the truncation error for polynomial approximation, 
p n (x), of a function /(x) is 

/(x) - P „(x) = ( x - f°Mf ~ ~ z ^ r+'((), («) 

{n + 1 j! 

where £ lies between the smallest and the largest x,. If one wants to minimize 
the error due to the term 


(x - x 0 )(x - Xj)...(x — x„) 


then the n+ 1 sample points should be chosen as the zeros of Chebyshev poly- 
nomial T n+ i(x), see [7]. This selection of grid points ensures that the error 
has the equal-ripple property which is characteristic of Chebyshev polynomi- 
als. A Chebyshev spectral method interpolates the highest order polynomial 
possible onto the n + 1 degrees-of-freedom defined by the point values of a 
function at the n + 1 zeros of T n+ i(x). The question to be addressed now 
is what if the grid is defined as the zeros of T n +i(x) but the polynomial is 
of lower order. That is, n could be, say, 128 whereas the polynomial on 
this grid could be of order 16. It is well-known that high order polynomial 
interpolation on uniform grids is essentially an ill-posed problem. 


4 Wavelet-based Grid and Order Selection 

The previous section introduced the idea of building very high order algebraically- 
generated difference operators on Chebyshev grids as a way of obtaining very 
high accuracy which is almost spectral in nature. This section will explore 
the idea of performing wavelet-based grid refinement on these Chebyshev 
grids as a way to obtain the necessary Chebyshev grid distribution near a 
boundary while having the ability to refine the grid away from the boundary 
for proper physical-space function resolution. 
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4.1 A Short Review of Wavelets 

To define Daubechies-based wavelets, see [6] for the original work and see 
[19] for an introduction to wavelet-based signal processing, consider the two 
functions 4>(x), the scaling function, and ip(x), the wavelet. The scaling 
function is the solution of the dilation equation, 

L - 1 

<t>{x) = V2 £ M( 2* - *), (47) 

k= 0 

where <f>(x) is normalized <f>(x)dx = 1, and the wavelet t/)(x) is defined in 
terms of the scaling function, 

L - 1 

xl>(x) = y/2^2 9 k<j>{^x - k). (48) 

k-0 

One builds an orthonormal basis from <j>(x) and ip(x) by dilating and 
translating to get the following functions: 

ti(x) = 2-i<f>(2~ j x-k), (49) 

and 

rpl(x) = 2~^ip(2~ j x - fc), (50) 

where j, k € Z. j is the dilation parameter and k is the translation pa- 
rameter. The coefficients H = {h k }^Zl and G = {gk)kZo are related by 
9k — (— 1 ) k hL-.k for k = 0, ..., L — 1. All wavelet properties are specified 
through the parameters H and G. If one’s data is defined on a continuous 

domain such as f(x) where x 6 R is a real number then one uses <f>{(x) and 

ti’k(x) to perform the wavelet analysis. If, on the other hand, one’s data is 
defined on a discrete domain such as f(i) where i G Z is an integer then the 
data is analyzed, or filtered, with the coefficients H and G. In either case, 
the scaling function <f>(x) and its defining coefficients H detect localized low 
frequency information, i.e., they are low-pass filters (LPF), and the wavelet 
4>{x) and its defining coefficients G detect localized high frequency informa- 
tion, i.e., they are high-pass filters (HPF). Specifically, H and G are chosen 
so that dilations and translations of the wavelet, iftKx), form an orthonormal 
basis of L 2 (R) and so that ip(x) has M vanishing moments which determines 
the accuracy. In other words, xj} 3 k { x ) will satisfy 

hi fijm = [ ip k (x)il>™(x)dx, (51) 

J —oo 
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where Ski is the Kronecker delta function, and the accuracy is specified by 
requiring that tp(x) — ^q(x) satisfy 



x m dx = 0, 


(52) 


for m = 0, M — 1 . Under the conditions of the previous two equations, for 
any function f(x) € L 2 (R) there exists a set {<£,-*} such that 


/(*) = d jk^K x ), ( 53 ) 

j£Z k£Z 


where 


/ OO 

f(x)ij>i(x)dx. 

-OO 


(54) 


The two sets of coefficients H and G are known as quadrature mirror 
filters. For Daubechies wavelets the number of coefficients in H and G, or 
the length of the filters H and G , denoted by Z,, is related to the number of 
vanishing moments M by 2 M = L. For example, the famous Haar wavelet 
is found by defining H as ho = h\ = 1. For this filter, H, the solution to 
the dilation equation (47), <f>{x ), is the box function: <j>{x) = 1 for x £ [0, 1] 
and <f>{x) = 0 otherwise. The Haar function is very useful as a learning 
tool, but because of its low order of approximation accuracy and lack of 
differentiability it is of limited use as a basis set. The coefficients H needed 
to define compactly supported wavelets with a higher degree of regularity 
can be found in [6]. As is expected, the regularity increases with the support 
of the wavelet. The usual notation to denote a Daubechies-based wavelet 
defined by coefficients H of length L is Dl- 

It is usual to let the spaces spanned by ^(x) and t/>£(x) over the parameter 
k, with j fixed, be denoted by Vj and Wj respectively, 


^ = span ii(z ), 
3 kez n 

(55) 

w > = T" 

(56) 

The spaces V’, and Wj are related by, 


... C Vi C Vo C V-i C 

(57) 
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and 


^ = V )+1 ©WO +1 , (58) 

where the notation Vo = Vi 0 W\ indicates that the vectors in Vi are orthog- 
onal to the vectors in Wj and the space Vo is simply decomposed into these 
two component subspaces. 

The previously stated condition that the wavelets form an orthonormal 
basis of L 2 (R) can now be written as, 

L\R) = ®W i . (59) 

Two final properties of the spaces Vj are that, 

nv,- {o}, (60) 

and 

U Vi = L 2 {R). (61) 

jez 

4.2 Grid Refinement on Uniform Grids 

The idea of using wavelets to generate numerical grids began with the ob- 
servation in [12] that the essence of an adaptive wavelet-Galerkin method 
is nothing more than a finite difference method with grid refinement. So, 
instead of letting the magnitude of wavelet coefficients choose which basis 
functions to use in a Galerkin approach, let the same coefficients choose 
which grid points to use and then think of the wavelet method in a colloca- 
tion sense. 

In other words, suppose a calculation begins with N evenly-spaced sam- 
ples of a function / and that some quadrature method produces N scaling 
function coefficients on the finest scale denoted by Vo- If the spacing between 
adjacent values in the vector / is Ax then this is also the physical-space 
resolution of any calculation done in Vo- Now, decompose Vo once to get 
Vo = Vi ® W\. Similarly speaking, the physical space resolution of Vi is 
2Ax and the refinement from the 2Ax physical-space resolution to the Ax 
physical-space resolution is dictated by the wavelet coefficients in W\. This 
is the reasoning which led to WOFD and to the following subroutine which is 
at the heart of WOFD. The remainder of the paper is concerned with giving 
the reader an idea of how the WOFD grid refinement software works. 
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4.3 Grid Refinement on Chebyshev Grids 

Chebyshev grids are not evenly-spaced in physical space, but are evenly- 
spaced in angle. That is, a Chebyshev grid comes from Xj = cos(0j), j = 

0 , ...,7V, where the angle 9j = ^ is evenly-spaced. The above described 
refinement mechanism can now be applied to the uniform angle grid point 
values to define a new numerical grid. That is, all the above grid refinement 
machinery can be applied to Chebyshev grids where each subspace Vj will 
coincide with a uniform angle or usual Chebyshev grid and each refinement 
subspace Wj will coincide with additional points being added to the usual 
Chebyshev grid. It is well known that the Chebyshev grid is the best grid, in 
terms of minimal error, for algebraic polynomial interpolation. A refinement, 
W\, on Chebyshev grid, Vi, to get V 0 = V\ ® W\ is disigned to begin with 
a grid which in one sense in perfect, the Chebyshev grid, and perturb from 
this grid. 

5 A New Numerical Method: WOFD2 

In [15] a numerical method was defined which was called the Wavelet- Optimized 
Finite Difference method or WOFD. WOFD used wavelets in their finite 
difference form. In essence, this meant that wavelets were used to choose 
a numerical grid and all computations were performed on this grid using 
arbitrary-grid finite difference operators. 

WOFD2 is an extension of this idea. WOFD2 uses wavelets to choose not 
only a numerical grid but also the order of the difference operator used on this 
grid. In addition, W0FD2 uses very high order finite difference operators 
on the order of 8, 16 or maybe 32. Furthermore, the physical-space grids 
are no longer evenly-spaced at every resolution but are Chebyshev. That is, 
wavelet-based grid generation, see [16], requires that a grid be selected from 
a uniform finest grid. But, high order polynomials can be highly-oscillatory 
on uniform grids. Therefore, WOFD2 works with Chebyshev grids at each 
resolution level. Recall, that Chebyshev grids x, = cos(0,-) are not uniform 
in the physical space variable x, but are uniform in the angle variable 6,. It 
is in this uniform angle variable 6, that grid refinement is performed. 

Using the grid selection mechanism in [16] to select order is a minor exten- 
sion of the idea that wavelets are very good at finding regions of the domain 
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at which a large numerical error is likely to occur. Numerical error will deter- 
mined by the truncation error of a polynomial which is locally interpolated 
to the data. The truncation error will be the product of intervals and a con- 
stant. Imagine the intervals are all a multiple of a smallest interval Ax then 
the key component in the truncation error will be (Ax) n . This component 
can be decreased by either decreasing the size of Ax or by increasing the 
order of the scheme, i.e., increasing n. Or, one can decrease Ax and increase 
n simultaneously. 

5.1 Comparison with hp- Refinement 

In the finite element literature, see [20], [1], [11], the idea of re fini ng the grid 
and increasing the polynomial order is known as hp-refinement. The theory 
from hp-refinement can certainly be applied to the new method WOFD2, 
even though WOFD2 works with polynomials only for generation of finite 
difference operators to be applied in the physical space. One of the most 
important results from hp-refinement theory from which WOFD2 can ben- 
efit is that when the function being differentiated is smooth then the rate 
of convergence is controlled by the polynomial degree. For the purpose of 
pulse propagation in aeroacoustics it is apparent that a high order differen- 
tiation, i.e., high order polynomial interpolation, will propagate the pulse 
more faithfully than grid refinement on the same pulse, assuming the pulse 
is smooth. 

5.2 Numerical Experiments with WOFD2 

This section will provide the results from numerous numerical experiments 
performed with WOFD2. For all the numerical experiments in this section of 
the paper a Gaussian pulse enters the domain from the right-hand side and 
travels to the left. The governing equation is the 1 dimensional hyperbolic 
wave equation, 

U t (x,t) = U x (x,t), U(x, 0) = e -c ( I-ir ) 2 , (62) 

for some constant c. See Figure 1 for a plot of this initial condition. Note that 
the domain extends from 0 to n. The final time for all simulation is 7t/2. The 
simulation is stopped at this value because at 7t/2 the Chebyshev grid has a 
maximum spacing between grid points and hence a minimum resolution. 
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Grid 

Pts 

Order 
of Acc 

l 2 

Error 

LqO 

Error 

Final 

Time 

64 

8 

3.82 * 10- 3 

1.78 * 10 -2 

tt/2 

64 

16 

2.95 * 10- 4 

1.30 *1(T 3 

7t/ 2 

64 

32 

1.74 *10“ 5 

7.10 * 10~ 5 

7t/2 

64 

48 

2.74 * 10- 6 

1.28 * lO" 5 

7t/2 


Table 6: W0FD2 of Accuracies 8, 16, 32, and 48 


5.2.1 No Adaptation 

First we consider tbe case of very high order finite differencing on a Cheby- 
shev grid. The grid size is kept fixed at 64 points, and the order is increased 
from 8 to 48. The errors decrease in a nice and uniform manner. No unusual 
numerical oscillations occur. 


5.2.2 Adapting Grid Only: Order 8 Spatial Differencing 

In this subsection the order of the spatial differencing is kept fixed at order 8. 
No results are found for threshold values of 10 _1 and 10 -2 . This is because it 
seems to be a characteristic of adaptive methods that a very rough threshold 
value can degrade the performance of the method. It is better to start with 
threshold values less than or equal to 10 -3 . The first row of the table is the 
worst possible performance where no refinement is done and grid is 64 points, 
and the last row of the table is best possible performance where the grid is 128 
points. Note that the software is constructed to work with both periodic and 
non-periodic boundary conditions, so that when the boundary conditions are 
non-periodic the possible number of points becomes 2^ + 1 which includes the 
right-hand boundary point. For periodic boundary conditions the number of 
grid points is 2 N since the right-hand boundary point is equal to the first 
point on the left-hand boundary. 
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Grid 

Pts 

t S 

Grid 

Order 
of Acc 


L<x> 

Error 

Thresh 

Final 

Time 

128/64 

65 

8 


H£Q£S9 


7r/2 

128/64 

77 

8 

21 

g££Ql!ld 


7t/2 

128/64 

79 

8 

■ Bl 

Bl 



128/64 

82 

8 

Eg£gg£5fl 


mam 


128/64 

84 

8 

(gigQiSfl 


HSB 


128 

129 

8 



0.0 

7t/2 


Table 7: W0FD2 Adaptive Grid, but Accuracy fixed at 8 


5.2.3 Adapting Grid Only: Order 16 Spatial Differencing 

Much of what was said for the 8th order table above can be said here. The 
second row of the following table shows how the performance can be slightly 
degraded for relatively large threshold values. In this case the degradation 
occurs at the threshold value of 10~ 3 . This is a minor point. Generally 
speaking, just start with a smaller threshold value. 

5.2.4 Adapting the Grid and Order 

This final table is the culmination of the paper and an example of W0FD2 
with all options in use. The grid is adjusted between a maximum density of 
128 and a minimum density of 64. The order of accuracy is adjusted between 
a maximum order of 16 and a minimum order of 8. The error converge in a 
nice manner toward the minimum error which occurs at the maximum grid 
density of 128 and the maximum order of accuracy of 16. 

The usual Chebyshev grid is evenly-spaced in angle 6i = iir/N for i = 
0, ..., N. In the physical space the grid distribution is x t - = cos (6,) which 
is shaped like a semi-circle. When one applies the wavelet grid adaptation 
to this evenly-spaced 8i then obtains in the physical space the distribution 
X{ = cos(#,) in the portion of the domain away from the pulse, and the twice- 
as-dense grid distribution x, = cos(jV/(2iV)) in the portion of the domain 
near the pulse. Note that the grid is the usual Chebyshev grid near the 
boundary. It is only safely away from the boundary that the grid density 
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Grid 

Pts 


*/ 

Grid 


Order 
of Acc 


L 2 

Error 


2.95 * IQ -4 


" OO 

Error 


Thresh 


Final 

Time 


128/64 


65 


16 


1.57 * 10 -3 


1.30 * 10 -3 

7.31 * 10" 3 


100.0 


7t/2 


128/64 


81 


16 


10 


-3 


2.34 * lO" 4 
2.43 * 10~ 5 


x/2 


128/64 


87 


16 


8.16*10 


-4 


10 


-4 


7t/2 


128/64 


89 


16 


8.09 * 10 


-5 


10 


-5 


2.36 * lO" 6 


9.20 * 10~ 6 


7r/2 


128/64 


87 


16 


1.09* 10~ 6 


10~ 6 

' W 7 ' 


n r/2 


128/64 


91 


16 


3.18 *10" 7 


1.05 *10" 7 


4.66 * 10- 7 
2.24 * 10" 7 


7r/2 


128/64 


93 


16 


10 


-8 


4.96 * in- 8 


tt/2 


128 


129 


16 


0.0 


7t/2 


Table 8: WOFD2 Adaptive Grid, but Accuracy fixed at 16 


Grid 

Density 

tf 

Grid 

Order 
of Acc 

l 2 

Error 

Loo 

Error 

Thresh 

Final 

Time 

64 

65 

8 

khhes 


100.0 

7t/2 

128/64 

81 

16/8 

Qj223Qj2l 


msm 


128/64 

80 

16/8 



10 -4 


128/64 

82 

16/8 


i 

mm i 


128/64 

84 

16/8 


ESSES! 

■ESI 


128/64 

86 

16/8 


EUllggi 



128/64 

87 

16/8 


i Bfl 

10" 8 


128/64 

90 

16/8 

5.12 * 10~ 8 


10" 9 


128 

129 

16 


2.24 * lO" 7 




Table 9: WOFD2 with Grid and Order Adaptation 


27 













































































makes an abrupt change in density. See Figure 2 for an example of an initial 
grid. 

If the numerical scheme is working properly then the pulse will propagate 
to the middle of the domain and be similar in shape to the initial condition. 
The best measure of this similarity is the L 0 Q error. At the final time the 
pulse will appear as in Figure 3. 

The Chebyshev grid is naturally more dense near the boundaries than in 
the middle of the domain. With the wavelet adaptation of this Chebyshev 
grid, the grid points can be kept dense while maintaining a the Chebyshev 
distribution throughout most of the domain. Again, the most important 
region of the domain for a Chebyshev distribution is near the boundary. See 
Figure 4 for the grid distribution at the final time when the pulse has reached 
the middle of the domain. 

Without grid refinement or order refinement the peak numerical error at 
the final time should be near the peak of the pulse, since it is this portion of 
the function which is most difficult to represent by polynomial interpolation. 
See Figure (5) for an example of such an error. 

If the wavelet refinement threshold is not sufficiently low then one will 
see the peak error appear near a region of the domain where there is a grid 
or stencil discontinuity. A ‘sufficiently low’ refinement threshold will on the 
order of the L <*, error when no grid or order refinement is executed. In 
Figure 6 noise is amplified at the interface where both the stencil and grid 
are refined. If the wavelet refinement threshold is adjusted to a smaller value 
then one can obtain an error profile similar to that in Figure (5). 

When both the stencil and grid are changed throughout the calculation, 
one finds a relatively wide stencil near the peak value of the pulse. For the 
example of a 17 point stencil with accuracy of 16 at the pulse and a 9 point, 
accuracy 8, stencil away from the pulse see Figure (7). 


6 Conclusion 

This paper has covered many topics related to the construction of a very 
high order adaptive order and adaptive grid numerical method which has 
been named the Wavelet- Optimized Difference Method 2, or W0FD2. First 
it was necessary to explore the various ways in which difference operators can 
be constructed. This included a comparison of difference operators generated 
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from algebraic, trigonometric, exponential, and cosine polynomials. Next, 
which type of polynomial would be best for the construction of very high 
order numerical differencing. The conclusion, which is not a big surprise, 
is that one should use algebraic polynomials on Chebyshev grids. The next 
step was to apply wavelet grid and order adaptation in order to be able 
to reduce errors throughout the domain by either increasing the order of 
the numerical method or by increasing the grid density in the appropriate 
region. The results of the numerical tests were very positive and it appears 
that WOFD2 will applicable to a large range of numerical problems. The 
version of W0FD2 which has been presented here has been ‘tweaked’ very 
little. That is, it worked essentially for the first time it was tried. This is 
encouraging because most high order numerical methods require some kind of 
filtering or other refinement. Future plans for W0FD2 would be, perhaps, to 
try to find a proof of stability and to apply the method in higher dimensions. 
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Figure 3: Pulse at Final Time 
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Figure 5: Typical Error at Pulse Peak 
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